R Notebook

This R notebook conducts various tests using the method outlined in Toward automated extraction and characterization of scaling regions in dynamical systems.

# Dependencies
library(ggplot2)
library(LaplacesDemon)
source("./generic.R")
source("./deshMethod.R")

states_visited <- c()
acfs <- c()

# Specific to autoregression 
a = 0.95
b = 2
exp_m <- 1 + 2*(a/(1-a))
steps = 10000

get_next_step <- function(prev) {
  w <- rnorm(n=1, mean=0, sd=1)
  next_val <- a * prev + b * w
  return(next_val)
}

get_init_state <- function( ) {
  init_sd <- b^2 / (1-a^2)
  init_state <- rnorm(n=1, mean=0, sd=init_sd)
  return(init_state)
}

# example method with one simulation using one set of p, q
# also plots histograms here. 
n <- 10 # chosen
p <- 3 # chosen
q <- 3 # chosen 
steps = 10000
a <- 0.9
r <- 0

states_visited <- asim(get_init_state(), get_next_step, steps)
acf_df <- c_est_series(states_visited)

# limiting the search to a finite number, since takes a long time computationally
stop_search <- which(acf_df$acf <= 0)[1]
stop_search <- 1.5 * stop_search

y <- get_log(acf_df$acf)[1:stop_search]
x <- c(0:(stop_search-1))
lim <- length(y)
results_ar1 <- desh_method(x,y, lim, n)

[1] 58
Time difference of 0.2918248 secs

new_weights <- get_desh_weights(results_ar1,  p, q, r)

[1] 1128

vars_ar1 <- get_slopes_sides(results_ar1, new_weights, TRUE)

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

vars_ar1

exp_m <- 1 + 2*(a/(1-a))
sprintf("Resulting RHS : %f, expected time: %f", vars_ar1$best_rhs, exp_m)

[1] "Resulting RHS : 20.932509, expected time: 19.000000"

This is the big test– does the method for many different a’s, p’s, and q’s.

# Generates simulations for each possible a we want to test, storing results in a dictionary. 
#poss_as <- c(-0.9, -0.7, -0.5, 0.5, 0.7, 0.9, 0.99)
poss_as <- c(0.5, 0.7, 0.9, 0.99)
steps = 10000
dict <- c()
for(a in poss_as) {
  states_visited <- asim(get_init_state(), get_next_step, steps)
  acf_df <- c_est_series(states_visited)
  dict[as.character(a)] <- acf_df
}

# Then we try the method on each simulation.
# Start by finding the autocorrelation functions. 
poss_ps <- c(1, 2, 3)
poss_qs <- c(1, 2, 3)
poss_rs <- c(0)

old_ms <- c()
new_ms <- c()
ps <- c()
qs <- c()
as <- c()
real_m <- c()
results_ar1s <- list()

for(a in poss_as) {
  acf <- dict[as.character(a)]
  acf_df <- data.frame(acf)
  colnames(acf_df) <- c("acf")
  stop_search <- which(acf_df$acf <= 0)[1]
  stop_search <- 2 * stop_search
      
  if(a < 0) {
    n <- 3
  } else if(stop_search > 500) {
    stop_search <- 500
  } else {
    n <- 10
  }
  
  if(stop_search < 5) {
    stop_search <- 10
  }
  y <- get_log(acf_df$acf)[1:stop_search]
  x <- c(0:(stop_search-1))
  lim <- length(y)
  results_ar1 <- desh_method(x,y, lim, n)
  results_ar1s <- append(results_ar1s, results_ar1)
}

[1] 20
Time difference of 0.0153501 secs
[1] 36
Time difference of 0.1190212 secs
[1] 160
Time difference of 6.823408 secs
[1] 500
Time difference of 8.277671 mins

# Then try the method on the given autocorrelation functions.
old_ms <- c()
new_ms <- c()
ps <- c()
qs <- c()
as <- c()
real_m <- c()
rs <- c()

for(a in poss_as) {
  print(a)
  start_time <- Sys.time()
  c <- 4
  M_upper_bound <- 6000
  triangle_df <- find_m(acf_df, c, M_upper_bound)
  m_suggest <- triangle_df[triangle_df$poss_Ms >= triangle_df$ct_ints,]
  m_suggest <- m_suggest$poss_Ms[1]
  m_suggest <- t_int(m_suggest, acf_df)
  
  acf <- dict[as.character(a)]
  acf_df <- data.frame(acf)
  colnames(acf_df) <- c("acf")
  stop_search <- which(acf_df$acf <= 0)[1]
  stop_search <- 1.5 * stop_search
      
  if(a < 0) {
    n <- 3
  } else if(stop_search > 500) {
    stop_search <- 500
  } #else {
    #n <- 10
  #}
  
  if(stop_search < 5 ) {
    stop_search <- 10
  }
  y <- get_log(acf_df$acf)[1:stop_search]
  x <- c(0:(stop_search-1))
  lim <- length(y)
  results_ar1 <- desh_method(x,y, lim, n)
  
  
  for(p in poss_ps) {
    for(q in poss_qs) {
       for(r in poss_rs) {
        exp_m <- 1 + 2*(a/(1-a))
        real_m <- append(real_m, exp_m)
        old_ms <- append(old_ms, m_suggest)

        new_weights <- get_desh_weights(results_ar1,  p, q, r)
        vars_ar1 <- get_slopes_sides(results_ar1, new_weights, plot_title=sprintf("a:%f, p:%f, q:%f", a, p, q)) #automatically plots RHS 
        new_ms <- append(new_ms, vars_ar1$best_rhs)
        
        ps <- append(ps, p)
        qs <- append(qs, q)
        as <- append(as, a)
        rs <- append(rs, r)
       }
    }
  }
  end_time <- Sys.time()
  print(end_time - start_time)
}

[1] 0.5
[1] 15
Time difference of 0.003687143 secs
[1] 10

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 10

[1] 10

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 10

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 10

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 10

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 10

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 10

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 10

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

Time difference of 59.63924 secs
[1] 0.7
[1] 27
Time difference of 0.04512501 secs
[1] 136

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

Time difference of 59.47326 secs
[1] 0.9
[1] 120
Time difference of 2.51026 secs
[1] 5995

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

Time difference of 1.037575 mins
[1] 0.99
[1] 500
Time difference of 7.548814 mins
[1] 119805

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805

Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

Time difference of 11.9294 mins

# Formatting and storing results
results <- data.frame(old_ms, new_ms, ps, qs, as, real_m)
colnames(results) <- c("Old Method", "New method", "p", "q", "a", "Real Time")
write.csv(results, "./results.csv", row.names = FALSE)
results

# plotting the densities for the best p, q for each a
for(a in poss_as) {
  results_a <- results[results$a == a, ]
  results_a$err_vs_real <- (results_a["Real Time"] - results_a["New method"]) / results_a["Real Time"]
  selected_vals <- results_a[which.min(abs(results_a$err_vs_real$`Real Time`)),]
  acf <- dict[as.character(a)]
  acf_df <- data.frame(acf)
  colnames(acf_df) <- c("acf")
  
  stop_search <- which(acf_df$acf <= 0)[1]
  stop_search <- 1.5 * stop_search
  
  if(a < 0) {
    n <- 3
  } else if(stop_search > 500) {
    stop_search <- 500
  } else {
    n <- 10
  }
  
  y <- get_log(acf_df$acf)[1:stop_search]
  x <- c(0:(stop_search-1))
  lim <- length(y)
  
  results_ar1 <- desh_method(x,y, lim, n)
  print("yo")
  print(results_ar1)
  new_weights <- get_desh_weights(results_ar1,  selected_vals$p, selected_vals$q)
  
  vars_ar1 <- get_slopes_sides(results_ar1, new_weights, plot=TRUE)
  
  
  
  
}

[1] 3
Time difference of 0.0002748966 secs
[1] "yo"

[1] 0

Error in density.default(results$slopes, bw = "nrd", adjust = 0.25, kernel = "gaussian",  : 
  argument 'x' must be numeric